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We calculate the efficiency of a rejection-free dynamic Monte Carlo method for d-dimensional 
off-lattice homogeneous particles interacting through a repulsive power-law potential r~^ . Theo- 
retically we find the algorithmic efficiency in the limit of low temperatures and/or high densities 

p+2 _d 

is asymptotically proportional to p 2 T 2 with the particle density p and the temperature T. 
Dynamic Monte Carlo simulations are performed in 1-, 2- and 3-dimensional systems with different 
powers p, and the results agree with the theoretical predictions. 
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I. INTRODUCTION 

Since their introduction in 1953 fl'| classical Monte 
Carlo (MC) methods mave matured into a useful tool 
for studying many different phenomena in different fields 
such as material science, high energy physics, and biology 
0) S 3 • For the study of the statics of a given model, the 
MC method can be viewed as a probabilistic method of 
performing multi-dimensional integrals 5] that could cor- 
respond to the partition (or grand partition) function of 
the model [^ . Various methods to increase the efficiency 
or accuracy of the method exist, including importance 
sampling 0, which allows the MC method to provide 
an estimate for the ratio of two integrals, for example to 
give an estimate of the average energy {E) of the system. 
Other advanced algorithms, such as Swendsen-Wang or 
cluster algorithms 13, iq, can also be used to alleviate 
difficulties associated with critical slowing down near T^ 
or being frozen into valleys for T < Tc. All of these 
advanced methods are allowed because they provide es- 
timates for the underlying integral(s) in an efficient fash- 
ion. For any of these MC methods, the system to be 
studied is in some configuration, which has a particular 
energy Ei, and an algorithm to obtain a new configura- 
tion j from configuration i is implemented. For M gener- 
ated configurations, the estimate for the average energy 
isgivenhy{E) = jjj:Z,E.,. 

If one is interested in the physical time evolution of 
a model system, the MC method can still be used. Al- 
though in principle any MC algorithm for statics could 
be used to study the time development through phase 
space of the model system, only certain MC methods 



will correspond to the actual time evolution of the sys- 
tem being modeled. In other words, many of the meth- 
ods mentioned above, such as Swendsen-Wang or cluster 
algorithms [2, Q, change the rate at which the system 
moves through phase space and consequently would usu- 
ally not be associated with the actual time development 
of the physical system. The older MC literature simply 
refers to the method as a MC method (see the refer- 
ences in Sec. 3.4 of Ref. [J]), even when the time devel- 
opment of the MC algorithm is assumed to correspond to 
that of the actual model system Q . More recently, use 
of MC methods to study the physical time dependence 
of a model system has been called either dynamic MC 
or kinetic MC. Although these two terms are sometimes 
used interchangably, there is an emerging distinction be- 
tween them. Kinetic MC has become the standard name 
for the case where physical time development is studied 
with known rate constants for the system to evolve from 
one state to another Q. These rate constants may be 
approximated under certain assumptions (such as appli- 
cability of transition state theory to atomistic systems) 
using ab initio methods ^ . (Use of rate constants also 
allows the transition from discrete MC steps to contin- 
uous time.) However, there are other instances where 
the physical time evolution of the system is desired while 
rate constants might be unavailable (for example per- 
haps transition state theory does not apply or an im- 
portant complicated multi-particle motion might be dif- 
ficult to conceptualize or calculate). In such cases the 
physical time development may sometimes still be de- 
rived from the underlying physical system, for example 
by studying the underl ying quantum mechanism for time 
development [3, [lO, [ill [I4I , or devising a method equiva- 



lent to the time development of the underlying equations 

MM- 

Such studies, which we concentrate on in this arti- 
cle, are called dynamic MC studies to distinguish them 
from static (equilibrium) MC or from kinetic MC studies. 
Thus we use the term dynamic MC in the same way as 
the recent book by Krauth [1^ . 

Frequently in dynamic simulations we need to work 
with long time scales at very low temperatures or in a 
strong external field. In these cases the standard dynamic 
MC method becomes very inefficient due to the high re- 
jection rate which requires a large number of trial moves 
before a change is made to the state of the system. Most 
advanced algorithms, such as Swendsen-Wang or clus- 
ter algorithms [2, , change the dynamic of the system 
thereby changing the time development of the system, 
which makes it impossible to study systems where the 
MC move is based on physical processes. The rejection- 
free MC (RFMC) method was proposed to overcome this 
problem with standard dynamic MC. The RFMC was 
first applied to discrete spin systems [1^1 ll3i llM| includ- 
ing the kinetic Ising Model [1]. It was later generalized 
to classical spin systems with continuous degrees of free- 
dom 19]. The RFMC method allows us to efhciently 
simulate a system with a high rejection rate without 
any changes of the original dynamics, since it shares the 
original Markov chain with the standard dynamic MC 
method. The RFMC method in this paper is very similar 
to the method labeled 'faster-than-the-clock algorithm' 
in Ref. M- 

The RFMC requires, to proceed by one algorithmic 
step, the values of all the probabilities of choosing a new 
state. Therefore, the computational cost of one step is 
larger than that of the standard MC. For this reason it is 
necessary to have a method to calculate its efficiency on 
a particular problem without implementing the method 
directly. Watanabe et al. developed a method [20| to cal- 
culate the efficiency of the RFMC for spin systems and 
hard particle systems. In the present paper, we evalu- 
ate the efficiency of the RFMC method for dynamic MC 
studies of d-dimensional particle systems with the par- 
ticles interacting through a repulsive short-range power 
law potential. Even though the bookkeeping involved in 
actually implementing the RFMC method may be sub- 
stantial, leading to more computer time per algorithmic 
step than the standard dynamic MC method, at temper- 
atures low enough or fields high enough where rejection 
rates are extremely high, the RFMC will be more efficient 
than the standard dynamic MC algorithm. 

The paper is organized as follows. In Sec. II, we present 
a review of the standard dynamic MC algorithm and the 
RFMC algorithm. In Sec. Ill, we provide analytical esti- 
mates for the efficiency of the RFMC method for repul- 
sive power law potentials. In Sec. IV, we show the results 
of our simulations in 1, 2 and 3 dimensions. Sec. V is de- 
voted to discussions and conclusions. 



II. DYNAMIC MONTE CARLO 

A. Standard Dynamic Monte Carlo 

The standard dynamic MC algorithm for particle sys- 
tems involves the following six iterative steps, with one 
iteration being called a MC step (MCS). We have used 
the term dynamic MC for this algorithm, rather than 
the term time-quantified MC used in I13. [l4| . since it 
has been used previously [lO, Ull [l3, ^M- I* is im- 
portant to remember that in dynamic MC the time in 
MCS is proportional to the physical time in seconds 
d, [10, Ell, [12, [11 [3 • The algorithm satisfies detailed 
balance. It is very similar to the time- qua ntification of 
the dynamic MC for Brownian ratchets [14| . 

1. Choose one particle randomly from the N particles, 
the chosen particle i is located at position roid,i. 

2. Choose a new position of the chosen particle ran- 
domly as Tncw.i = T'o\d,i + Ar, with Ar chosen uni- 
formly over a d-dimensional hyper-spherical volume 
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r( 



1 



(1) 



with a radius rchoose and the gamma function F. 
The probability density for choosing the new posi- 
tion of the chosen particle is df-Xi/Vchoosc- 

3. Reject the new position if it is located outside the 
'cage' formed by line segments joining its nearest- 
neighbor (nn) particles. (In the usual way, the 
'cage' is defined in this off-lattice simulation on the 
basis of a Voronoi diagram (or Delaunay triangu- 
lation), but can be often equally well defined for 
our homogeneous high-density systems as particles 
within a certain distance of the chosen atom 20].) 

4. Evaluate the energy difference AE'i = E'new.i — 
^-oid.i between the new and the old positions of the 
chosen particle i. 

5. Decide whether to accept the trial move by com- 
paring a random number with the move probability 
which is a function of AE".;. For example, we can 
use the Metropolis criteria to choose the transition 
probability as 



P{ 



^ncw,i I '^o\d,i) 



1 if AE, < 0, 

exp(— /3A£'i) otherwise. 



(2) 



6. If the trial move is accepted, move the particle to its 
new position, otherwise leave it in its old position. 



B. Dynamic Monte Carlo Without Rejections 

The rejection- free MC (RFMC) method was developed 
to overcome the decrease in the efficiency of the standard 



dynamic MC method in cases where the rejection rate is 
high, for example, in systems at a very low temperature. 
The efficiency of the standard dynamic MC algorithm 
is the rate at which it changes the current state of the 
system, this is the fraction of accepted moves to the total 
number of move attempts. One algorithmic step of the 
RFMC involves the following procedures. 

1. Compute the time to leave the current state (the 
waiting time iwait)- This is the number of trial 
states which would be rejected in the standard dy- 
namic MC. Hence in one algorithmic step the time 
is advanced by twait . 

2. Advance the time of the system by iwait- 

3. Calculate probabilities A^ for each of the N parti- 
cles, where Xi denotes the probability that the trial 
move of particle i would be rejected in the standard 
dynamic MC given that it was the particle chosen 
for the trail move. Explicitly 



A,: - 1 - 



1 



V, 



cho 
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■ro\d,i)d'^Xn 



(3) 



4. Then choose a particle with the probability pro- 
portional to 1 — Ai, that is, the probability that 

the particle i is chosen is given by — rr • 

Ef=i(i-A.) 

Therefore, the particle which is easy to move has a 
higher probability to be chosen. 

5. Choose a new position of the chosen particle. This 
is accomplished using the probability density 



P(r„ 



fo\d,i)d''-x„ 



V, 



choose 



(1 - A.) 



(4) 



The above procedure does not contain any rejection step, 
and therefore, each RFMC algorithmic step always in- 
volves a change of state of the system. 

The efficiency of the RFMC method is inversely pro- 
portional to the rejection rate of the standard dynamic 
MC, but they share the same dynamics. Therefore the 
efficiency of the RFMC is related to the inefficiency of 
the standard dynamic MC. For this reason the efficiency 



of the RFMC will be proportional to i„ait which is given 
by 0,[2ll: 



^wait — 



Inf 



1. 



(5) 



Here f is a random number uniformly distributed on 
(0, 1], [-J is the integer part, and A = -^ X]i=i Ai is the 
probability to stay in the current state (that the move 
to the trial state will be rejected) after one standard dy- 
namic MC trial move. The units of time are in MCS, but 
can be quantified with physical time |a.[l(|JllL[l3.[l3l.[l4]|. 
We have used the Metropolis method [l| as shown in 
Eq. (m in this example and in our simulations in the next 
section. However, the RFMC algorithm would also work 
for different functional probabilities such as the Glauber 
or heat-bath dynamic [2, y, l9| or a phonon dynamic 

MM. 



III. EFFICIENCY OF RFMC FOR POWER-LAW 
POTENTIALS 

Consider d-dimensional particles with a repulsive 
power-law potential and A^nn nearest-neighbors (nn) [the 
particles that form its 'cage']. The potential between any 
two nn particles is 



V{r) 



(7r-(^r^^ 

r> 



ro 



r >ro 



(6) 



where r is the distance between particles, p is the power, 
ro is the cut-off distance and a is the length that gives 
the strength of the interaction, respectively. This po- 
tential function represents the hard repulsive core of any 
potential, such as the Lennard-Jones potential which has 
p = 12, and this repulsive part becomes dominant at high 
densities. 

We have chosen the origin of the coordinate system 
to be at the position of the chosen atom i. The energy 
difference of the atom i is given by AEi = Ui (x) — Ui (0) 
with X the trial position. The efficiency of the rejection- 
free algorithm at low temperatures and/or high densities 
is given for the Metropolis dynamic by 



J 



{exp[-pAE]) 



r(f + i) 



'^^ ^choose "'-°° 



d'^xO 



cage 



exp • 



{-/3[f/,(f)-[/,(0)]}, 



(7) 



where the interaction energy for particle i is given by interatomic potential 
the power-law dependence, p, of the repulsive part of the 



Wn 



JV„ 



fe = l 



(8) 



with Xk the position of the k nn atom of the chosen since Vchooso — 2rchoosc, -'Vnn 
atom i. Here 0cagc restricts the integrand to be non 



2, and 



zero only with the cage formed by the nn particles. The 
angular brackets denote an average over all allowed states 
of the system weighted with the Boltzman weight at each 
configuration. Since we are interested in the system at 
high densities or at low temperatures, we can utilize the 
Laplace saddle-point integration approximation [22| 



d''xP{x)^P{xo)^[^ (9) 



where the integrand P{x) is strongly peaked around x 
xq and 
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A-.- = — 



-In [P{x)] 



(10) 



We assumed the chosen particle is at or near its local 
energy minimum, i.e., P (xq) w 1. Therefore 



{exp[-(3AE]) 



r(i + i)ri 



choose 



\ 



i2nr 



det A 



(11) 



since after making the derivation in Eq. (|10p we get a 
factor of /3 and the only values that are not are when 
i — j we have defined Aij = (3Aij, and the determinant 
of A is thus proportional to (3'^. This immediately gives 
that the temperature dependence of the average waiting 
time is proportional to (3^ — T^2 , 

Furthermore, because of the power-law approximation 
of Eq. ([8]) , and that two partial derivatives must be taken 
for the saddle-point approximation Aij ^ ''nn '^ where 
Tnn is the nn distance if all nn atoms are equidistant from 
the chosen atom. The particle density p is proportional 



' nn ' 



Equation ([7|), therefore becomes 
(exp [-PAE]) 



[r^r^] 



T^ 



choose r^ 



(12) 



and the average time between acceptances in the dynamic 
MC procedure is 



{U 



1 



choose r 

(exp[-/3A£;]) '^ Ti 



(13) 



Equation p3|) is the main result of this paper. The result 
is very general, both for various dimensional particles and 
for various power laws, as well as being general for the 
explicit dynamic that is used in the MC procedure. 

Throughout the present paper, we assume all the 
atoms have identical potentials. In the following, we 
calculate the explicit expression of Eq. P^ for several 
conditions. For the 1-dimensional system, the average 
waiting time is given by 



(i- 



wait/f^^l 



■^^choosc 

^(P+2)/2 
' nn 



aPp{p -1-1) (0 2 



T-K 



Vt' 



(14) 



U,{x)-US) 



dx^ 



2aPp{p+l) 

p+2 
''nn 



(15) 



For the 2-dimcnsional system which has a hexagonal lat- 
tice as the ground state, the average waiting time is given 
by 



(t 



wait/rf=2 



Sc^^'P^'-^hoosc P' 



2TrP+^ 



since Vchooso = Trr^, 



A,, 



choose' ^nn 



U,{x) ~ [/,(0) 



nn 

6, and 



dxidxi 



T ' 



3crV 

'J ,,P+2 ' 
' nn 



(16) 



(17) 



with the Kronecker delta Sij . For the 3-dimensional sys- 
tem which has the face-centered-cubic (FCC) lattice as 
the ground state, the average waiting time is given by 



,, V _ Kp(P-I)]^ ^choose 



twait/d^S 



P 



10+3P 3 - 
3y/TT2 4 i 2rn 



Ti ' 



(18) 



since Vchooso == (4/3)7rr3|^^^^^, A^nn = 12, and 



Aij — 



C/,(x)-[/,(0) 



dxidxi 



= -5r 



aPp{p- 1)2^ 



„P+2 
' nn 



^^0 



(19) 

Note that, the density and the temperature dependence 
in a simple-cubic lattice is equivalent to that of Eq. (|19p . 
while the coefficient is different since iVnn — 6 and An = 
-%(aMp-l))/< 



rP+2 
nn ■ 



IV. SIMULATIONS IN d = 1, 2, AND 3 

Following the methodology described above, we per- 
formed simulations for 1-, 2-, and 3-dimensional systems. 
The goal is to locate where the asymptotic results of 
Eq. ^ hold for our RFMC method. The density of the 
system p is defined to be p = N{2a/LY, with the num- 
ber of particles iV, the radius of the particles a, the linear 
system size L, and the dimensionality of the system d, re- 
spectively. Throughout all the simulations rchooso is set 
to 0.05 and the cut-off radius Tq is set to tq = l.lr„„(p), 
with r„„(p) the nearest- neighbor distance for the given 
density. The simulations were performed for four values 
of p: 2, 4, 6, and 12. For all values of d, the density is 
fixed at p = 2.0 in order to study the temperature de- 
pendence, the temperature is fixed at T = 0.001 to study 
the density dependence. All simulations were performed 
starting from the ground state, with periodic boundary 
conditions, with the simulated volume such that an inte- 
ger number of unit cells of the ground state fit into the 
volume. 
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FIG. 1: (Color online) (a) Temperature dependence of the 
average iwait in d = 1 for A^ = 200 particles, the solid lines 
are the power law fits, (b) Density dependence of fwait in 
d — 1, the solid lines are power law fits with p > 0.6. 



FIG. 2: (Color online) (a) Temperature dependence and (b) 
density dependence of the average tv,ait in d — 2 ioi N — 80 
particles. The solid lines are the power law fits. 



We study the 1-dimcnsional system with N = 200 
particles which are located on the line of length L with 
periodic boundary conditions. The temperature depen- 
dence of the average for t„ait is shown in Fig. [1] (a), 
in this case with statistics from 10^ MCS per particle 
(MCSp). The power law fit gives (t„ait) -- r-0-49(i) with 
the correlation coefficient r — 1.00(1) for p — 2,4 and 
6, while for p = 12 we obtain (twait) ^ T^^-^"'!) with 
r = 1.00(1). We have fit the region with T < 10"\ 
T < 10", T < 10\ and T < 10^ for respectively p = 2, 4, 
6, and 12. All of these results are in agreement with our 
prediction, i.e., (iwait) ~ T~"'^. The density dependence 
of iwait is shown in Fig. [T] (b) with the number of trials in 
this case 10^ MCSp, or lO'' MCSp for high density and 
p = 12. The power law fit in this case, all for p > 0.6, 
for p = 2 gives (t^ait) -^ p^-^^^^^ with r = -1.00(1), for 
p = 4 gives (iwait) -^ p2.99(i) ^j^j^ J. ^ -1.00(1), for 
p = 6 gives (iwait) -- p4.oo(i) ^j^-j^ ^ ^ -1.00(1), and for 
p = 12 gives (iwait) ~ p''-°°(i) with r = -1.00(1). Again 
our results agree with our asymptotic predictions, i.e., 

(twait) ^ p(P+2)/^ 

We study the two-dimensional system with A^ = 80 
particles distributed on a triangular lattice of length L 
and width \/3L/2 with periodic boundary conditions for 
both axes. Fig. [5] (a) shows the temperature dependence 
of the average iwait , in this case with the number of sam- 



ples 10^ MCSp. The power law fits, all for p > 1, give 
(t^ait) -- T-0-99(i) with r = 1.0(1), for p = 2, 4, 6 and 12. 
The density dependence is shown in Fig. [21(b), 10^ MCSp 
are taken for p = 2, 4, 6 and 10^ for p = 12 at high den- 
sities. The power law fit for p = 2 gives {t^ 



-1.00(1), for p = 4 gives (i^ait) 



^1.99(1) 
.3.00(1) 



p400(l) with 



with r = 

with r = -1.00(1), for p = 6 gives (iwait) 
r = -1.00(1), and for p = 12 gives (iwait^ ^ P 
with r — —1.00(1). These results show excellent agree- 
ment with the asymptotic prediction (twait 
Eq. C 



„7.09(1) 
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We also study the 3-dimensional system with A^ — 500 
particles located on an FCC lattice with periodic bound- 
ary conditions for all directions. The temperature de- 
pendence and density dependence of the average twait 
are shown in Fig. [3] (a) and Fig. [3] (b), respectively. 
The power law fits of the temperature dependence gives 
(twait) - T-i-49(i) with r = 1.00(1), for p = 2,4,6 and 
12. The power law fit, all for p > 1.1, for the den- 
sity dependence for p = 2 gives (twait) ^ 
r = -1.00(1), for p = 4 gives (twait) - 
r = -1.00(1), for p = 6 gives (twait) - 
r = -1.00(1), and for p = 12 gives (twait) 
r = —1.00(1). The results agree with the asymptotic 
predictions (twait) -- T-3/2p(p+2)/2. 



pl-98(l) with 

^3.02(1) ^itij 

p4-01(l) with 

. p6.96(l) .^^j^jj 





FIG. 3: (Color online) (a) Temperature dependence and (b) 
density dependence of the average fwait for a 3-dimensional 
FCC system with A^ — 500 particles. The solid lines are the 
power law fits. 



V. SUMMARY AND DISCUSSION 

We studied the efRciency of the Rejection Free Monte 
Carlo (RFMC) method for systems having particles inter- 
acting through repulsive power-law potentials. The den- 
sity and temperature dependence of the average waiting 
time has been predicted to be, 



{K 



p+2 
P 2 
d 
T2 



(20) 



with the dimensionality of the system d, density p, and 
the temperature T, respectively. These theoretical re- 
sults are valid asymptotically for large p and/or low T. 
Monte Carlo simulations were performed and the results 
showed good agreement with the asymptotic prediction. 



This study shows how efScient the RFMC method is in 
low temperature or in high density regimes. Assume 
the wall-clock time per algorithmic step for the standard 
dynamic Monte Carlo algorithm is t^, and the average 
wall-clock time per RFMC algorithmic step is ti (both 
of which are expected to be almost independent of T and 
p). Because of the extra bookkeeping involved in pro- 
gramming the RFMC method, ii > to- Nevertheless, the 
RFMC method will be more efficient (use less wall-clock 
time) whenever ii < io(^wait)- This inequality will always 
be satisfied for low enough T or high enough density. 

The RFMC method does not change the dynamic of 
the MC move associated with the underlying physical 
dynamics, and therefore makes possible the study of sys- 
tems with a fixed physical dynamic. It is very important 
keeping the dynamics unchanged, since the change in the 
dynamics of the MC move can cause a str ong influence 
in certain dynamic physical properties [ll|, ll3, ^M, 1231 • 

Although we studied a repulsive core power-law po- 
tential, we expect the equivalent behavior of the waiting 
time for more realistic potentials, as long as we work with 
the system at low temperatures and/or high densities. 
A further avenue of study, therefore, could be to calcu- 
late the efficiency of the RFMC method for more general 
potentials such as Lennard- Jones, or those derived from 
density functional theory. A related study could also be 
of the efficiency of the RFMC in the case of two or more 
types of particles. The actual implementation and uti- 
lization of the RFMC method in particle simulations can 
be attempted now that the ultimate behavior of the al- 
gorithmic efficiency has been determined. 
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